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Background: To extract information about the neutrino properties from the study of neutrinoless 
double-beta {QvfifS) decay one needs a precise computation of the nuclear matrix elements (NMEs) 
associated with this process. Approaches based on the Shell Model (ShM) are among the nuclear 
structure methods used for their computation. ShM better incorporates the nucleon correlations, 
but have to face the problem of the large model spaces and computational resources. 
Purpose: The goal is to develop a new, fast algorithm and the associated computing code for effi- 
cient calculation of the two-body matrix elements (TBMEs) of the Qu/5/3 decay transition operator, 
which are necessary to calculate the NMEs. This would allow us to extend the ShM calculations for 
double-beta decays to larger model spaces, of about 9-10 major harmonic oscillator shells. 
Method: The improvement of our code consists in a faster calculation of the radial matrix ele- 
ments. Their computation normally requires the numerical evaluation of two-dimensional integrals: 
one over the coordinate space and the other over the momentum space. By rearranging the ex- 
pressions of the radial matrix elements, the integration over the coordinate space can be performed 
analytically, thus the computation reduces to sum up a small number of integrals over momentum. 
Results: Our results for the NMEs are in a good agreement with similar results from literature, 
while we find a significant reduction of the computation time for TBMEs, by a factor of about 30, 
as compared with our previous code that uses two-dimensional integrals. 

Conclusions: We developed a new, fast, efficient code for computing the TBMEs that are used to 
calculate the NMEs necessary for the analysis of the 0^/3/3 decays. A rearrangement of the expres- 
sions of the radial matrix elements allow us to perform only one integration instead of computing 
two-dimensional integrals. This leads to a significant reduction of the computational time, which 
makes us confident that it is now possible to rapidly, accurately, and efficiently calculate the TBMEs 
for many major harmonic oscillator shells. 

PACS numbers: 23.40.Bw, 21.60.Cs, 23.40.-s, 14.60. Pq 
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I. INTRODUCTION 



The neutrinoless double beta (OV/3/3) decay is a be- 
yond Standard Model (SM) process of major interest for 
understanding the neutrino properties. Indeed, its dis- 
covery would decide if neutrinos are their own antipar- 
ticles PQ, and would give a hint on the scale of their 
absolute masses. That is why there are intensive inves- 
tigations on this process, both theoretical and experi- 
mental. The present status of these investigations can 
be found in several more recent reviews [2]-[S], which 
also contain therein a comprehensive list of references 
in the domain. Of particular interest is the effective neu- 
trino mass, a parameter entering the OV/3/3 decay half- 
lives, which depends on the neutrino masses, neutrino 
oscillating parameters and Majorana phases. Until now 
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this decay mode has not yet been confirmed by indepen- 
dent measurements and thus, one can only extract upper 
limits of this parameter from the existent experimental 
lower limits of the half-lives. However, to do this we 
also need a precise computation of the nuclear matrix el- 
ements (NMEs) which also enter the half-lives formula. 
An accurate calculation of these NMEs is one of the most 
important challenges in the theoretical study of the 0^/3/3 
decay. 

Typical calculations of the NMEs are performed us- 
ing a bare transition operator [5]. This is almost always 
the case even if one uses different approaches: pnQRPA 
0-dD], Shell model(ShM)[12]-[16], IBA [mill], PHFB 
[It?] and energy density functional (EDS) method [2U] - 
[2T] , which are the most common methods of calculation 
of these matrix elements. This is equally true even if 
one uses an improved transition operator that considers 
higher order effects in the nucleon current (HOC) [TT1I22) . 
In principle the most reliable of these approaches to per- 
form calculations for the NMEs (relevant for OV/3/3 de- 
cay) is the ShM, since it incorporates all types of corre- 
lations and uses effective nucleon-nucleon (NN) interac- 
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tions which are checked with other spectroscopic calcula- 
tions for nuclei from the same region. However, it has to 
face the problem of the large model spaces and the associ- 
ated computational resources. Also, it is well known that 
in ShM calculations of the two-neutrino (2^/3/3) matrix el- 
ements the Gamow- Teller operator needs to be quenched, 
to better describe the experimental data for beta decays 
and charge-exchange reactions. Therefore, it is impor- 
tant to know if the Or//3/3 transition operator has to be 
effectively modified when used in relatively small model 
spaces. Work in this direction was recently reported in 
Ref. [53] where an effective operator was analyzed for the 
0^/3/3 decay of 82 Se in the jj44 model space consisting of 
the /g/2) £>3/2i Pi/2; 59/2 orbitals. For these calculations 
up to 8 major harmonic oscillator shells (MHOS) were 
used, which implies that one needs all two body matrix 
elements of the 0^/3/3 transition operator in these large 
spaces. In addition, there were recent proposals [24J [25] 
to investigate the modifications of the transition operator 
in increasingly larger shell model spaces for a fictitious 
0^/3/3 decay of a p— shell nucleus. 

The calculations reported in Ref. [53] were performed 
using a bare operator without higher order contributions 
in the nucleon current. In these calculations the inte- 
gral over momentum in the transition operator can be 
analytically done, which makes the calculation of its two- 
body matrix elements very fast. It is however known that 
the effect of the higher order contribution in the nucleon 
current is a reduction of the Ov/3/3 matrix element by 20- 
30 %. This reduction could be further amplified by the 
equivalent effective operator. Therefore, it is important 
to investigate this effect, which would require knowledge 
of the two-body matrix elements of the bare transition 
operator in a large number of MHOS, e.g. 8 to 12. 

In our previous works [T5], [2S], we started to develop 
an efficient nuclear ShM approach to accurately calculate 
the NMEs for both 2v/3(3 and 0^/3/3 decay modes. The 
approach used in Ref. [TB] to calculate the two-body ma- 
trix elements (TBME) of the transition operator that in- 
cludes higher order terms in the nucleon current needs to 
calculate two-dimensional integrals, on the relative mo- 
mentum and the relative coordinate. This approach was 
sufficiently fast for calculating the two-body matrix ele- 
ments in a single major shell, such as pf -shell. However, 
calculations of these two-body matrix elements in 8-12 
major shells would be intractable with this approach. 

In this paper we present a new improved (fast, effi- 
cient) ShM code which reduces substantially the comput- 
ing time of calculation of the two-body matrix elements 
of the transition operators for the 0^/3/3 decay. In a sim- 
pler version of the code when finite nucleon size (FNS) 
and higher order nucleon currents correction (HOC) ef- 
fects are not included, the integrals over the neutrino po- 
tentials can be performed only in the coordinate space, 
as has been done in Ref. [12]. In the full version of the 
code, where these effects and the short range correlations 
(SRC) effects are included, normally two-dimensional in- 
tegrations need to be done, one in the coordinate space 



and one in momentum space. The main improvement 
in this code is a rearrangement of the two-body matrix 
elements that allows us to do the radial integrals (the in- 
tegrals in coordinate space) analytically when harmonic 
oscillator (HO) single particle wave functions are used. 
Therefore, only the integration over the momentum re- 
mains to be performed numerically. We first compare 
our results for NMEs with other similar results from lit- 
erature performed with both ShM and other methods. 
Then, we compare the CPU times of our code with the 
CPU times of our previous code [TB] , for the same calcu- 
lations. We note that these times decrease significantly. 
We get an estimation of an average CPU time per TBME 
and note that the new code proves very promising for 
more elaborate calculations in many MHOS. The paper 
is organized as follows: in section 2 we describe our for- 
malism giving the relevant expressions for the radial in- 
tegrals. In section 3 we present and discuss our results 
and the last section is devoted to the conclusions. We 
complete the paper with an Appendix where expressions 
of several quantities used in the formalism described in 
section 2 are given. 



II. DESCRIPTION OF THE FORMALISM 

The Ov/3/3 decay (Z, A) ->■ (Z+2, A)+2e~ requires that 
the neutrino and the antincutrino are identical and mas- 
sive particles. Considering that this decay occurs only 
by exchange of light neutrinos between nucleons and in 
the presence of left-handed weak interactions, the lifetime 
can be expressed as: 

(r^ 2 y l = G°»(E ,Z) I M°» | 2 , (1) 

G 01 ' is the phase space factor depending on the energy 
decay Eq and nuclear charge Z, and (m„) is the effec- 
tive neutrino mass parameter depending on the first row 
elements of the neutrino mixing matrix U e {, Majorana 
phases e lai and the absolute neutrino mass eigenstates 
rrii (see e.g. Ref. [5]). The nuclear matrix elements are: 

M°» = M^ T ■ Mf , (2) 

where M% y T and Mf are the Gamow-Teller (GT) and the 
Fermi(F') parts, respectively. Usually a tensor part ap- 
pears as well, but the numerical calculations have shown 
that its contribution is small [16] : consequently, it will 
be neglected in the following. The GT and F parts are 
defined as follows: 

i\C = £(0+||T_ m r_ n O^||0+) , (3) 

where 0" m are transition operators (a = GT, F) and the 
summation is over all the nucleon states. 
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Due to the two-body nature of the transition opera- 
tor, the matrix elements are reduced to sum of products 
of two-body transition densities (TBTD) and matrix el- 
ements for two-particle states (TBME), 



M 0v = 



]T TBTD {j p j p , 

(jpjp'', Jtt 1 1 T- IT- 2 II jnjn' J Jtt) , (4) 



The two-body transition operators Of 2 can be expressed 
in a factorizcd form as: 



0? 2 = N a sM ■ i& fc > 



(5) 



where N a is a numerical factor including the coupling 
constants, and and are operators acting on the 
spin and relative wave functions of two-particle states. 
Thus, the calculation of the matrix elements of these 
operators can be decomposed into products of reduced 
matrix elements within the two subspaces |16j . The ex- 
pressions of the two-body transition operators are: 



n GT 
u 12 



a\ ■ a 2 H(r) 



H(r) 



(6) 



The most difficult is the computation of the radial part 
of the two-body transition operators, which contains the 
neutrino potential. We will refer to its computation in 
more detail. Neutrino potential is of Coulomb type, de- 
pending weakly on the intermediate states, and is defined 
by integrals of momentum carried by the virtual neutrino 
exchanged between the two nucleons pj] 



H a (r) = 



2R 



3G 



Jo(qr) 

w 1 

1 

jo(qr)V a (q)q 2 dq 



1 



(E) 



q 2 dq 



(7) 



where R = 1.2A 1 / 3 fin, u = \J q 2 + to 2 is the neutrino 
energy and jo(<Z r ) is the spherical Bessel function. We 
use the closure approximation in our calculations, and 
(E) represents the average excitation energy of the states 
in the intermediate odd-odd nucleus, that contribute to 
the decay. The expressions of h a (a = F, GT) are 



(8) 



and 



h GT (q 2 



G\{q 2 ) 

2G 2 M (« 2 
3 



3q 2 



Ami 



(9) 



where m T is the pion mass, m p is the proton mass and 

G M (<7 2 ) = (^-Mn)Gy( 9 2 ), (10) 
with ([j,p — /!„) = 4.71. 



The expression ([9| includes finite nucleon size (FNS) 
and higher order terms in the nucleon currents (HOC). 
The Gy and Ga form factors which takes into account 
the FNS effects are: 



= 9v 



A 2 



A 2 , + q 7 



(11) 

For the vector and axial coupling constants we used gy = 
1 and gA — 1.25, and the values for the vector and axial 
vectors form factors we used are Ay = 850MeV and 
A A = 1086MeV 0. 

For computing the radial matrix elements 
(nl\H a \n' V) we use the harmonic oscillator HO wave 
functions ip n i(lr) and ip n i p (r) corrected by a factor 
[1 + f(r)], which takes into account the short range cor- 
relations induced by the nuclear interaction: 

ipni(r)^[l + f(r)}iP„i(r) . (12) 

For the correlation function we take the functional form 



f(r) 



-c ■ e 



l-br 2 



(13) 



where a, b and c are constants which have particular 
values for in different parameterizations |30j-|31|. as it 
will be discussed in the next section. 

As an observation, we note that in the limit m v = — > 
uj = q and when HOC and FNS effects are neglected, 
the integral over q in Eq. (7) can be easily done and the 
expression of the neutrino potential becomes 

H a cx - [sin«J5) r) Ci((E) r) - cosUE) r) Si((E) r)l (14) 
r 

where Si(z) and Ci(z) are the sine and cosine integral 
functions [TS]- Then, the calculation of the radial inte- 
grals over the neutrino potential reduces to a single inte- 
gral in the coordinate space, which can give a hint about 
the magnitude of the NMES within an error of about 25- 
35%. This approximation was used in Ref. |12j where 
the NMEs for 48 Ca are calculated. 

Including HOC and FNS effects the radial matrix ele- 
ments of the neutrino potentials become: 

/>oo 

(nl | H a (r) | n'V) = / r 2 dr^ nl (rty„'i'(r) [1 + /(r)] 2 
Jo 



x / q 2 dqV a (q)j (qr) 
Jo 



(15) 



where v is the oscillator constant. 

As one can see, if all the nuclear effects are included, 



the calculation of the radial integrals (151 requires the 



numerical computation of two integrals, one over the co- 
ordinate space and the other over the momentum space. 
However, one can reduce the computation to only one 
integral proceeding as follows: one can rearrange the ex- 
pression of the radial integral in coordinate space as a 
sum of terms with the same power of r. The dependence 
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of r appears from the product of the HO wave functions, 
the correlation function and Bessel functions. First, one 
can write the product of two HO wave functions as a 
sum over the terms with the same power in r: 



n+n 



~>Pni(r)ipn'i'{r) = ^ At +l ' +2s (nl,n'l') 



x (2 I ,) i±i ^ ± V w y +r+2s , (16) 

where j4; + ;< +2 s ar e coefficients independent of r whose 
expressions are given in Appendix. Then, one adds the 
contribution of the factor [1 + f(r)] 2 that brings depen- 
dence of r in powers of 0, 2 and 4: 



[l + /(r)f = l-2ttf 



2bce~ ar V + c 2 e- 2ar 



- 2bc 2 e- 2ar r 2 + b 2 c 2 e- 2ar V (17) 

Finally, the computation of the radial matrix elements 
requires to perform integrals of the form: 



/>oo 

X Q (/x;m) = / q 2 dqV a (q) 
Jo 









^ r 


x 




) 2 (2*0 










Jo 



dr e-» r r m j (qr) (18) 



where = j/, j/ + a, v + 2a and m is integer. In this 
expression the integration over r can be done analytically 
and one gets: 



(21/)" 



¥-1 



x(m-l)!!^(-l) 

fc=0 

Thus, Z a (fi;m) becomes: 



k I 2 



- 1 
k 



e 

(2fc + l)!!(2 / Lt) fc9 



,2A' 



(19) 



J Q (/u;m) 



fe=0 



fc 2 



Jafak) (20) 



where J^ a {fi; k) are integrals over momentum: 

Ja(M; k) = 



{2k + 1)!! (2/i) fc 



(21) 



Finally, the radial matrix element can be expressed as a 
sum of integrals over the momentum space: 

n+n 

(nl | H a (r) | n'l') = ^ A I+r+a .(nf,n / Z / )X: a (m) (22) 

s=0 

where JC a (m) is a sum of six I Q (/x; m) integrals over mo- 
mentum. Its expression is given in Appendix. 



III. NUMERICAL RESULTS AND 
DISCUSSIONS 

We developed a new code for computing the TBMEs 
necessary for the ShM calculations of the NME involved 
in Of/3/3 decays, based on the formalism described in sec- 
tion 2 that reduces the two-dimensional integrals in Eq 
to sums of one-dimensional integrals of Eq. (21 1 



The new code has a flexible user interface, which allows 
the selection of various nuclear effects. As it was shown 
in the previous section, the improvement consists in a 
simplification in the computation of the radial matrix el- 
ements which allows us to perform the integration only in 
the momentum space. We first check our code comparing 
our results for the total NMEs, for 48 Ca and 82 Se with 
similar results from Refs. [12] and [8]. In order to obtain 
the M 0l/ nuclear matrix element, the two-body transi- 
tion densities are calculated using the method described 
in Ref. [IS]. For 48 Ca we used GXPF1A [27] effective 
interaction in the full pf model space, and for 82 Se we 
used JUN-45 [5S] effective interactions in the jj'44 model 
space. 

As one can see from Table I, our results are in good 
agreement with previous ones, provided that the same 
nuclear nuclear effects are included in the calculations. 
Our code enables us to include in a flexible manner the 
nuclear effects, such as SRC of Jastrow [^3] type with 
Miller-Spencer and coupled cluster model (CCM) [30] 
with Argonne V18 and CD-Bonn [30]-[31] parameteri- 
zations, FNS, and HOC. The SRC parameters entering 
Eq. (13) are the same as in Table II of Ref. [T5]. For 



comparison, we used in our calculation the same nuclear 
effects as those used in the corresponding references. 





48 C7a 


* 2 Se 


^ present work 


0.573 


2.47 


[E] (2010 ShM) 


0.57 




[13] (2008 ISM) 


0.59 


2.11 


P3] (2009 ISM) 


0.61 


2.30 


[S] (2007 QRPA) 




2.77 



TABLE I: Comparison between the results of the present work 
(*' and other similar results from the references indicated. In 
the calculation we used SRC of Jastrow type, FNS and HOC. 



We have also analyzed the performance of our code in 
getting an improved computing speed. In Figure [l] we 
show the single-core CPU times needed to compute the 
TBMEs. 

In the case of 48 Ca, there are 94 TBMEs requiring a 
total of 6.29s of CPU-time on our test machine equipped 
with Intel Xeon X5560 CPUs. This translates into an 
average of 6.7 • 10 _2 s for each individual TBME. When 
computing the product of wave functions, the dependence 
on the n and I quantum numbers of the nucleon orbits is 
reflected in the CPU-times, as one can see in the differ- 
ence between the average CPU-time/TBME of 48 Ca and 
those of 82 Se. 82 Se has required a total of 6.30s for the 
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Single-threaded CPU-times on an Intel Xeon CPU X5560 @> 2.80GHz 




8 Ca 82 Se 
FIG. 1: CPU-times for the computation of the TBMEs 



computation of its 65 TBMEs, thus needing an average 
time of 9.7 • 10~ 2 s for each TBME. Even then, we can 
still calculate TBMEs for 82 Se almost as fast as for the 
simpler case of 48 Ca nucleus. 

Figure [l] also shows the contribution of the SRC and 
FNS+HOC effects to the total computation time for the 
TBMEs. "Bare" means that neither SRC nor FNS + 
HOC effects were considered. With our new method and 
code, we obtain an improvement in speed by a factor of 
about 30, as compared to the code used in Ref. [16], 
where more than three minutes where needed instead of 
6.3 seconds. 

The performance of the new code makes us confident 
that it is now possible to rapidly, accurately and effi- 
ciently compute TBMEs for many nuclear shells. This 
task is very challenging for the TBME code of Ref. [TB] . 
For example, if one wants to investigate the effective tran- 
sition operator in only 8 MHOS [23] , one needs to calcu- 
late about 434k TBMEs (GT plus F). The actual average 
time/TBME is about 1.7 seconds, but as remarked above, 
is increasing with the raise of the angular momenta of the 
single particle orbits involved. Using a conservative es- 
timate of about 10 seconds/TBME one could conclude 
that one needs about 50 days of single-threaded process- 
ing power to calculate all necessary TBMEs. This time 
could be reduced by a factor of say 500 if the calculation 
of the TBMEs is distributed via a load-balancing algo- 
rithm [32 , when using 1000 cores with 50% efficiency. 
However, this reduction might not be sufficient if 9 or 10 
MHOS need to be used. The new algorithm presented 
here could be extremely useful in reducing the calcula- 
tion time by another factor of about 30. 



IV. CONCLUSIONS 

We developed a fast, efficient code for computing the 
TBMEs, which are part of the NMEs necessary for the 
analysis of the 0^/3/3 decays. The improvement consists 
in a faster computation of the radial matrix elements 
using correlated HO wave functions. Their computa- 



tion normally requires the numerical evaluation of two- 
dimensional integrals, one over the coordinate space and 
the other over the momentum space. By rearranging the 
expressions of the radial matrix elements, the radial inte- 
grals can be performed analytically over the coordinate 
space, thus the computation reduces to sum up a small 
number of integrals over momentum. We check our code 
by comparing the values of the NMEs for 48 C 'a and 82 Se 
calculated with our new code with similar results from lit- 
erature and we found a quite good agreement. Further, 
we estimated the CPU-times for one single core needed 
to compute the TBMEs with our code and compare them 
with the similar CPU-times obtained with our previous 
code requiring two-dimensional integrals. We find a sig- 
nificant reduction of the computational time, by a actor 
of about 30. We also estimated the average CPU-time per 
single TBME in the cases 48 Ca and 82 Se and found very 
small values. This achievement makes us confident that it 
is now possible to rapidly accurately and efficiently com- 
pute TBMEs for many major harmonic oscillator shells, 
which were very time-consuming in our earlier approach. 
The calculation of the TBMEs in 8 MHOS could be done 
in about 1-2 days using the present single-threaded code. 
Extension to more than 8 MHOS would require the par- 
allelization of the code using a load-balancing algorithm. 
These TBMEs can be further used to investigate the ef- 
fective transition operator needed for Oi//3/3 decay analy- 
ses. 



V. APPENDIX 



Al. The HO radial wave functions are given by: 



ipnl(r) = N n i exp 



r l L^Kr 2 , (23) 



where v is the oscillator constant, N n i is the normaliza- 
tion constant 



N„ 



2 n nl 



(2Z + 2ra + l)!! 



(2i/) 



2i + 3 / 2 



(24) 



(i+ L ) 

and Ln 2 (vr 2 ) is the Laguerre associated polynomials: 



Ln [UT > ~ 2^! 
x E (fc) ( 2 f + 2 1 fc + l)!! ^ 



(25) 



A2. The expression for the Ai + ii + 2s(nl, nil') used in Eq. 



G 



(221 is: 



Ai + i> +2a {nl,n'l') 



n\(2l + 2n + 1)!! n'\(2l' + 2n' + 1)!! 



2"' 



i 



kl(n- fe)!(2/ + 2fe + l)!! 



1 



fc / !(ra'-A; / )!(2Z / + 2fc / + l)!! 



(26) 



where a, 6, and c are the SRC parameters entering Eq. 
pl. 



with macc(0, s — n') < k < min(n, s) , A: + fc' = s 



A3. The expression for the IC a {m) used in Eq. (22) is 



JC a (m) = — \L a (y\ to) — 2cl a (v + a; to) 
V2j/ 

+ 2c j Z Q (^ + a;m + 2) + (?X a (y + 2a; to) 



(27) 



2c 2 ( — \ l a {v + 2a;m + 2) 
2^ 



<- | — | J a (i/ + 2a;m + 4)] , 
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